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We consider propagation models that describe the spreading of an attribute, called "damage", 
through the nodes of a random network. In some systems, the average fraction of nodes that remain 
undamaged vanishes in the large system limit, a phenomenon we refer to as exhaustive percolation. 
We derive scaling law exponents and exact results for the distribution of the number of undamaged 
nodes, valid for a broad class of random networks at the exhaustive percolation transition and in 
the exhaustive percolation regime. This class includes processes that determine the set of frozen 
nodes in random Boolean networks with indegree distributions that decay sufficiently rapidly with 
the number of inputs. Connections between our calculational methods and previous studies of 
percolation beginning from a single initial node are also pointed out. Central to our approach is 
the observation that key aspects of damage spreading on a random network are fully characterized 
by a single function specifying the probability that a given node will be damaged as a function of 
the fraction of damaged nodes. In addition to our analytical investigations of random networks, we 
present a numerical example of exhaustive percolation on a directed lattice. 

PACS numbers; 89.75.Da, 02.50.Ey, 02.10.Ox, 05.50.+q 



I. INTRODUCTION 
A. Overview 

Propagation models on lattices or more general graphs 
describe the spreading of some discrete signal through 
a set of discrete entities. In the most general terms, 
the signal corresponds to some qualitative change that 
causes the entity to interact differently with its neigh- 
bors. Examples include the spreading of damage in power 
grids [3, 0, the spreading of disease through a population 
I2I I^JnlL spreading of a computer virus on the Inter- 
net la, 13 1 or the alteration ofgene expression patterns in 
a cell due to a mutation [H In the general case, the 
individual entities are represented as nodes in a graph 
where the links indicate paths a long which the signal can 
spread [13 [II [11 [11 111 El [13^ Because the signal 
can be thought of as disrupting the static or dynamical 
state of the original system, we refer to its propagation 
as spreading damage, though in many cases the "dam- 
age" may enhance a desired property or simply represent 
some natural dynamical process. A single instance of a 
given spreading process initiated from a particular subset 
of nodes is often called an avalanche. 

In analyzing spreading processes, one is often inter- 
ested in the transition between those that die out quickly 
and those that spread to a finite fraction of the system 
in the large-system limit, a transition that may occur 
as the probability of transmitting damage across links 
is varied. This percolation transition is relevant for sys- 
tems in which the fraction of initially damaged nodes 
tends to zero in the limit of infinite system size. The 
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order parameter for the transition is the average fraction 
of nodes damaged in a single avalanche, which remains 
zero for small transmission probabilities and continuously 
increases when the probability rises above a threshold 
value. We will refer to this as the sparse percolation 
(sp) transition. The SP transition occurs for spreading 
processes in which the probability that a node becomes 
damaged is zero unless at least one of its neighbors is 
damaged. (If this probability were nonzero, a nonzero 
fraction of the nodes would always get damaged.) 

For a certain class of propagation models, there is an- 
other transition of interest. When the fraction of ini- 
tially affected nodes remains fixed as the system size is 
increased, the fraction of nodes that remain undamaged 
can undergo a transition from finite values to zero at 
transmission probabilities above some threshold. We re- 
fer to this as the exhaustive percolation (ep) transition. 
The EP transition occurs only for propagation models in 
which the probability of a node remaining undamaged is 
zero when all of its neighbors are damaged (all of its in- 
puts in the case of a directed graph) . We assume also that 
there is a nonzero probability for a node to remain un- 
damaged if it has at least one undamaged input. There is 
then one more condition for the ep transition: the density 
of directed loops of any specified size must vanish in the 
large system limit. For any loop there is a finite probabil- 
ity that no member of the loop will be damaged, since no 
member of the loop can have all of its inputs damaged 
until one of the members becomes damaged through a 
probabilistic event. Thus EP is not observable on spatial 
lattices of the type generally encountered in statistical 
mechanics, ep is observable, however, on directed lat- 
tices and on graphs in which the nodes serving as inputs 
to a given node are selected at random. 

In this paper we derive the probability distribution for 
the number of undamaged nodes at the EP transition on 
random graphs for a general class of propagation mod- 
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els exhibiting what we call unordered binary avalanches 
(uba). This is analogous to finding the distribution of 
avalanche sizes at the usual percolation transition, but 
here we are asking for the distribution of the number of 
nodes not participating in the avalanche. 

As an application of our EP results, we consider 
the problem of identifying unfrozen nodes in a random 
Boolean network (rbn). In a rbn, each node has a bi- 
nary state that is updated according to a rule that takes 
the values of some other nodes as inputs. The dynam- 
ics of RBNs has been investigated extensively; see, e.g., 
[11 111 111 113 in mm mill, a rbn can have sev- 
eral dynamical attractors, but some nodes might have 
the same value at all times on all attractors. Such nodes 
are called stable and the set of stable nodes is important 
for the dynamics in RBNs HE IH • 

Almost all stable nodes in a broad class of rbns can 
be identified through a dynamic process that was intro- 
duced by Flyvbjerg 26] and formalized to facilitate nu- 
meric simulations by Bilke and Sjunnesson |13|. We call 
the stable nodes that can be identified by this dynamic 
process frozen (and nodes that are not frozen are called 
unfrozen). Provided that the Boolean rule distribution is 
symmetric with respect to inversion of any subset of in- 
puts, the set of frozen nodes can be identified through an 
UBA in which frozen inputs cause new nodes to become 
frozen (damaged). Most rule distributions that have been 
examined in the literature exhibit this symmetry. The 
requirement is satisfied, for example, for any model that 
assigns given probability p for obtaining a 1 in each entry 
of the truth table for each node. 

This paper is organized as follows. We first develop 
the notation and basic definitions required for discussing 
UBAs in general. In Section llBI we give an introduction 
to the UBA formalism from the perspective of percolation 
processes. A more formal description is given in Sec- 
tion ^ followed by a numerical illustration of the basic 
concepts. In Section lllll we present analytic derivations 
for UBA in random networks with emphasis on ep and 
the EP transition. We also present explicit results for 
the special case of Erdos-Renyi networks with a natural 
choice for the avalanche rules. 

In Section llVl we show how to apply the UBA formalism 
to obtain the statistics of frozen nodes in two-input rbns. 
In the present context, this serves as an illustration of the 
general theory, but this particular example was also the 
primary motivation for studying EP. The results on rbns 
arc consistent with those found by Kaufman, Mihaljev, 
and Drossel. [13 . The main advantage of using the ep 
formalism for this problem is that it makes clear how 
the calculation can be extended to networks with more 
than two inputs per node, including networks with an in- 
degree distribution that (with a low probability) allows 
arbitrarily large in-degrees. 



B. Basic definitions 

An unordered binary avalanche (uba) is defined as a 
spreading process with the following properties: 

Binary states: the state of each node can be charac- 
terized as a binary variable s, with s = meaning 
undamaged and s = 1 meaning damaged; 

Boolean rules: the state of each node is determined by 
a Boolean function of the states of its input nodes; 

Order independence: the probability of having a 
given set of nodes damaged at the end of the pro- 
cess does not depend upon the order in which nodes 
are chosen for updating. 

Order independence refers to the dynamics of the 
spreading process or a simulation of it. In such a simula- 
tion, one typically chooses a site and updates it according 
to a rule depending on the states of sites that provide in- 
puts to it, repeating the process until a test of every site 
yields no change in the state of the system. We are inter- 
ested in cases where the order in which sites are chosen 
for possible updating has no bearing on the final state of 
the system. 

UBA is a natural extension of site or bond percolation. 
To determine the avalanche size distribution in site per- 
colation, for example, one identifies an initial subset of 
damaged sites and then tests neighbors of damaged sites 
to see whether the damage spreads to them. After a given 
site is tested for the first time, its value is permanently 
fixed. The process is iterated until no new damaged sites 
are generated. See, e.g., Ref. [s^]. This method of in- 
vestigating site percolation is equivalent to assigning all 
sites a value, then beginning with a damaged site and 
determining all of the damaged sites in a connected clus- 
ter. Site percolation where each site has the probability 
p to be occupied can be recast as a UBA system as fol- 
lows. Let each site be associated with a rule that is an 
OR-rule of all of its neighbors with probability p and is 
a constant with probability 1 — p. Then the above de- 
scribed site percolation is achieved by first selecting the 
rules and clamping the value of a given site to 1, and then 
repeatedly updating the system according to the Boolean 
rules. In this situation, the Is in the final state mark a 
site percolation cluster. A more practical way of simu- 
lating the same UBA is to determine probabilistically the 
Boolean rule at each site only when that site is first en- 
countered in the percolation process and to update only 
those nodes where the rules have been determined. 

To ensure order independence in UBA, it is sufficient 
to require that each Boolean function is non-decreasing, 
meaning that if one of the inputs to the rule changes from 
to 1, the output is not allowed to change from 1 to 0. 
For non-decreasing Boolean functions, if a specific node 
is eventually going to be assigned the value 1 during an 
avalanche, updating other nodes to 1 first cannot change 
the outcome. 
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We are particularly interested in UBAs that are initi- 
ated by damage at a set of nodes comprising a nonzero 
fraction of the total number of nodes. Such a process 
would be relevant, for example, if the probability that 
any given node is damaged at the start is independent of 
the system size. 



To clarify both the distinction between ep (exhaustive 
percolation) and SP (sparse percolation) and the similar- 
ities between them, we describe a particular case of a 
propagation model that exhibits both transitions. Con- 
sider a graph with a total of N nodes, some of which have 
three input links each while the others have no input links 
at all. The graph is random in that the node supplying 
the input value on any given link is selected at random, 
but stays fixed throughout the avalanche. Let vq be the 
fraction of nodes with no inputs. Define a spreading pro- 
cess as follows: The initial condition is that all nodes with 
no inputs are considered damaged. Each other node is 
now selected in turn to see whether the damage spreads 
to it. If a node has one damaged input, the probability 
that it will be damaged is pi ; if it has two damaged in- 
puts, the probability of damage is p2 (with p2 > Pi)', and 
nodes with three damaged inputs are guaranteed to be- 
come damaged (pa = 1). These probabilities are realized, 
for example, by the following Boolean rule distribution: 
a 3-input OR-rule with probability pi ; a 3-input majority 
rule with probability Pi; and a 3-input AND-rule with 
probability 1 — p2- 



As N goes to infinity, the number of initially dam- 
aged nodes can be a nonzero number that grows slower 
than N, meaning that fo goes to zero as N goes to in- 
finity. In this limit, the SP transition occurs at pi = 1/3 
and the spreading from each initially damaged node is 
described by a Galton-Watson process. In a Galton- 
Watson process, a tree is created by adding branches to 
existing nodes, with the number of branches emerging 
from each node drawn from a fixed probability distribu- 
tion. Such branching processes have been investigated 
extensively. (See, e.g., Rcf. 31].) In particular, the cor- 
respondence to Galton-Watson processes means that for 
critical SP, the probability of finding n damaged nodes 
scales like n"^/^ for K n < TV [l [s^l ■ 



For any nonzero value of j/q, the ep transition occurs 
forp2 satisfying (1— P2)(l~t'o) = 1/3 (assuming that this 
value of p2 is greater than pi.) The analysis described in 
Section lllTl Drovides a method of calculating the probabil- 
ity P{u) of having u undamaged nodes in this case. The 
result in the large N limit is P{u) ^ P{0)u~^/'^ for large 
u. A difference between EP and SP is that both P(0) and 
the cutoff on the distribution scale with for ep, 

while for SP only the cutoff scales with N . 



II. INTRODUCTION TO EXHAUSTIVE 
PERCOLATION 

A. Formal description of UBA 

We now describe a formalism and establish some no- 
tation that is suitable for a detailed treatment of UBA. 
Let TV denote the number of nodes in a network with 
a specified set of links and let the nodes be indexed by 
j = 1, . . . , A^. The network state is described by the vec- 
tor s = {si, . . . , Sat}. Let Kj denote the number of inputs 
to node j, and let kj denote the vector of Kj inputs to 
node j. Furthermore, let R denote a Boolean function 
and let Iij{R) denote the initial probability that node j 
has the rule R. [It is required that R has precisely Kj 
inputs for Hj{R) to be nonzero.] 

To efficiently simulate UBA, we keep track of the infor- 
mation that is known about each node at each step in 
the process. In particular, it is important to keep track 
of whether or not the change from to 1 of a given input 
has already been accounted for in determining the out- 
put. The simplest way to do this is to introduce an extra 
state 0* that labels a site whose rule R implies an output 
value of 1 but for which the update to 1 has not yet been 
implemented. When a node changes its state from to 
0*, it is a silent change in the sense that the Boolean 
rules at the other nodes treat an input 0* exactly the 
same as 0. To retrieve the final state of the network, all 
occurrences of 0* must be updated to 1. When a single 
update to 1 is made, the information that the given node 
has value 1 is passed along to all nodes with inputs from 
it. The values of these nodes may then change from to 
0*. The conditional probability that the value of node i 
is updated from to 0* when j changes value from 0* to 
1, is given by 



where is the value of k^ after Sj has been updated and 
Pi(ki) is the probability that Ril^i) — 1: 

Pi(kO = 5]i?(k,)n,(i?). (2) 

R 

The numerator in Eq. is the probability that Ri pro- 
duces a 1 after the update of node j minus the probability 
that Ri produced a 1 before the update. The denomina- 
tor is the probability that node i had the value before 
the update. 

Let 11^(1) denote the probability that the rule at node 
i has output 1 regardless of its input values. If some 
particular nodes are selected for initiation of the UBA, 
1141) is set to one for these nodes [which means Iii{R) — 
for all other rules]. 

We are now ready to present a formal algorithm for 
determining the final state of an instance of UBA on a 
finite network. We carry out the following procedure 
(where := denotes the assignment operator): 
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1. Sj := for all j; 

2. Sj := 0* with probability 11^(1) for each j; 

3. Some j with Sj = 0* is selected; 

4. Si := 0* with probability Ui{s,j) for each i with 

= 0; 

5. Sj 1; 

6. StepsOHSlare iterated as long as there exists a node 
in state 0*. 

UBA can also be considered on infinite networks, but 
that requires a more technical description of the process. 
First, the choices of j in step 3 for both descriptions must 
be such that any given j that satisfies the conditions in 
step 13 will be selected in a finite number of iterations. 
Second, the ensemble of final states needs to be defined 
in terms of a suitable limit process because the stopping 
criterion in step can not be applied to an infinite sys- 
tem. 

Note that the dynamics is only dependent on the prob- 
ability functions {Pi(ki)}. That is, the precise rule 
distributions affect the avalanche results only through 
their contributions to Pi . Because the Boolean rules are 
non-decreasing functions, Pi(ki) is also a non-decreasing 
function. In fact, every non-decreasing function, /(k^), 
with values in the interval [0, 1] can be realized by Pi(ki) 
for a suitable Boolean rule distribution. One such rule 
distribution can be constructed as follow: for each i and 
each ki, select a random number y from a uniform dis- 
tribution on the unit interval and set Ri(\ii) = 1 if and 
only if y < /(kj. 



B. An example of EP on a lattice 

To illustrate the concepts of UBA and ep, consider a di- 
rected network on a two-dimensional square lattice with 
periodic boundary conditions. Each node in the lattice 
has integral coordinates where i + j is odd and 

the node at (i, j) receives inputs from the two nodes at 
(i— 1, j±l). The rule for propagation of damage to a node 
is either OR or AND, with probabilities ^(i.j) (or) = r and 
(and) = 1 — r, respectively. 

Figure ^ displays an avalanche that is initiated by 
letting each node be initially damaged with probability 
p = 1/8. A node assigned OR becomes damaged if ei- 
ther of its neighbors one layer above is damaged; a node 
assigned and becomes damaged if and only if both neigh- 
bors above it are damaged. This means that 



^i(k(...)) = 



if k, 
if k, 
if k, 



(0,0) 



{(0,1), (1,0)} (3) 



(1,1) 



Note that clusters of damaged nodes formed in an 
avalanche initiated by a single damaged node cannot con- 
tain any holes, as the uppermost undamaged node in the 
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FIG. 1: An example of UBA on a lattice, displaying undam- 
aged nodes (dots), initially damaged nodes (filled circles), and 
nodes damaged during the avalanche (empty circles). Each 
node has either an OR-rule or an AND-rule with inputs from 
its neighbors in the row immediately above the node. The 
probability for a node to be initially damaged is p = 1/8 and 
the probability for obtaining an OR-rule is r = 0.3. Periodic 
boundary conditions are used and the first row and column 
are repeated in gray after the last row and column to illustrate 
the periodic boundary conditions. 



hole would have to have two damaged inputs and hence 
would become damaged when updated. 

For localized initial damage, the SP threshold is found 
at r = 1/2. Above this value of r, domains of dam- 
age tend to widen as the avalanche proceeds. Since the 
growing cluster has no holes, this is simultaneously an ep 
transition. The ep transition can be found for smaller 
values of r in lattices where each node is initially dam- 
aged with a given nonzero probability p. [For every ini- 
tially damaged node, Il^i j)(l) is set to 1, meaning that 
Pi(k(jj-)) = 1 for every value of kj-j j).] 

Figure|21shows the average number of unaffected nodes 
as a function of r for p = 1/8 on lattices with periodic 
boundary conditions. The numerics displayed in Fig. [21 
clearly suggest that there is a second-order ep phase tran- 
sition. Furthermore, these numerical results suggest that 
the avalanche in Fig.^is within the parameter regime for 
EP and that ep does not occur in this case due to finite 
size effects. 

For the case r = 0, it is possible to map the EP transi- 
tion onto ordinary, directed, site percolation on the same 
lattice. When all nodes in the lattice have AND-rules, the 
following algorithm may be used to determine whether a 
given node will be damaged: select a node; put a mark 
on the selected node unless it is initially damaged; and 
recursively mark each initially undamaged node that has 
an output to a marked node. The selected node will 
get damaged if and only if this recursion ends in a fi- 
nite number of steps. The algorithm describes ordinary 
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FIG. 2: The average fraction of undamaged nodes for UBA 
on a lattice of the type shown in Fig. as a function of 
the selection probability p for OR-rules and the probability 
p = 1/8 for initial damage. The lattice has periodic boundary 
conditions and covers a square that has a side of 10, 10^, 10^, 
and 10* lattice points, respectively, with steeper curves for 
larger systems. The statistical uncertainty in the estimated 
mean is less than the line width. 



directed site percolation where the initiaUy undamaged 
(damaged) nodes are considered active (inactive) sites 
and the process propagates in the opposite direction rel- 
ative to the UBA. We therefore expect the ep transition to 
occur for a value of p equal to 1 — Pc, where Pc = 0.70549 
is the threshold for directed site percolation |^ and we 
have confirmed this with numerical tests. Further study 
of EP on the lattice is beyond the scope of this paper. 



C. Suppression of EP by resistant motifs 

In the lattice example above, the fact that the network 
had no feedback loops smaller than the lattice size was 
important. In general, ep is suppressed by the presence 
of short feedback loops. As already noted, for EP to oc- 
cur, it is required that the output of each rule in the 
rule distribution is 1 if all of its inputs have the value 1. 
Otherwise, there would be a finite fraction of nodes that 
keep the value regardless of the influence from the rest 
of the network. Generalization of this reasoning allows 
us to rule out ep in other situations, indicating that EP 
is most likely to occur in directed or highly disordered 
networks. To pursue this idea, we introduce the notion 
of resistant motifs. 

A motif is a small network with a particular arrange- 
ment of internal links. A given motif may occur many 
times in a network with different rules assigned to its 
nodes and with different configurations of external in- 
puts. A motif is resistant with respect to a given en- 
semble of rule assignments if the probability of damage 



entering the motif when all external inputs are damaged 
is strictly less than unity. For the rule distributions that 
we consider for the ep transition in random networks, 
each node has a nonzero probability of being assigned a 
rule that sets its output to if at least one of its inputs 
is 0. Thus when all of the nodes in a feedback loop of 
any length have the value 0, there is a nonzero probabil- 
ity that they will all remain even if all external inputs 
to the loop are set to 1. Every feedback loop of a given 
length is therefore a resistant motif. 

If the number of occurrences of a resistant motif grows 
linearly with the network size, there will in total be a 
finite fraction of nodes that remain unaffected with a 
finite probability. For such networks, ep cannot occur 
in the limit of large systems. Examples include typically 
studied regular lattices and small world networks with 
link directions assigned so that short feedback loops are 
prevalent. 

The problem resistant motifs can be avoided in ran- 
dom networks having a mean indegree (K) that is well- 
defined and independent of N, in which case the number 
of feedback loops of a given length approaches a con- 
stant. Though the total number of resistant motifs may 
grow with system size, the larger motifs have a low prob- 
ability of avoiding damage. For large N, the out-degree 
distribution is a Poisson distribution with a mean value 
of (K). The outputs emerging from a given node form 
a tree with approximately (isT)'" nodes at the mth level. 
Thus, the probability for a given node to be part of a cy- 
cle of m nodes is approximately {K)™/N, which means 
that the typical number of feedback loops of length m 
is approximately {K)"^ /m. On the other hand, the loop 
may contain either initially damaged nodes or some rules 
that allow damage to enter from external inputs. The 
probability that this will not occur decays exponentially 
with m. If the decay is faster than (if)"™, the density of 
nodes in undamaged resistant motifs will approach zero. 

In summary, EP (for the considered type of rule dis- 
tributions) is excluded on lattices with a high density of 
feedback loops. For random networks, however, the frac- 
tion of nodes in undamaged resistant motifs can go to 
zero in the large N limit. This property allows ep to oc- 
cur on random networks as demonstrated in the following 
section. 



III. EP ON RANDOM NETWORKS 
A. Criteria for EP 

Consider a network such that the inputs to each node 
are chosen randomly and uniformly from all nodes in the 
network and the probability functions {Pi(ki)} are de- 
termined from a given distribution of Boolean rules. For 
such networks, UBA can be handled analytically. 

Define g{x) as the probability for a rule in the ran- 
dom network to output 1 if each input has the value 1 
with probability x. The function g reflects the probabil- 
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ity for propagation of damage to a single node, for the 
considered instance of UBA. We refer to g as the dam- 
age propagation function. In random networks, Pi(ki) 
is independent of i and can be replaced by Pi(k). Let 
K denote the number of components of k, i.e., the num- 
ber of inputs to the considered node. g{x) can then be 
expressed as 



9{x) 



E 

K=0 



ke{o,i}^^ 



a;^(l-a;)^-^Pi(k), (4) 



where I is the number of Is in k and P{K) is the prob- 
ability to draw a rule with K inputs. 

Let N denote the total number of nodes, and let no, 
no* , and ni denote the number of nodes with the values 
0, 0*, and 1, respectively. With these definitions and 
the fact that Ui{s,j) is independent of i for the random 
network, the role of {Pi(ki)} is taken over by g{ni/N) 
and Eq. is replaced by 



g{n'jN) ~ gjnjN) 
l~g{ni/N) ' 



(5) 



where n[ — ni 



1 . This means that the size of the net- 
work, the number of initially damaged nodes, and the 
damage propagation function g taken together are suffi- 
cient to uniquely determine the stochastic spreading pro- 
cess. 

After one pass of the update steps (from Sec- 
tion III All , the new values Uq and n^n. of no and no* are 
given by 



no = no 



and 



where 



nn* 



no* 



(5-1 
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(6) 



(7) 



(8) 



with Bn{a) being a stochastic function that returns the 
number of selected items among n items if the selection 
probability for each of them is a. The avalanche ends 
when no* — 0. 

The number of damaged nodes, n, in a complete 
avalanche is the final value of ni, whereas the number 
of undamaged nodes, u, is the final value of no. An order 
parameter for the system is (f) = ImiN-tooin/N), where 
the average is taken over the ensemble of networks. The 
SP transition is found when (j) changes from zero to a 
nonzero value, whereas the EP transition is found when 
(j) reaches 1. 

To understand the typical development of an 
avalanche, it is convenient to change from the variables 
no, no*, and ni, which are constrained to sum to N, to 
the variables xi = ni/N and 



no 



l-5(a;i)' 



(9) 



As long as no* > 0, the average value of c after a single 
update is given by 



1 - 
no 



5W) 



1 - .gW) 



(10) 

(11) 

(12) 



Hence, as long as no* > for all members of an ensemble 
of avalanches, (c) (the average of c over the ensemble) is 
conserved as the avalanche proceeds. 

From Eqs. Q-® and the definition of c, the variance 
in c can be calculated. We begin by computing the in- 
crement of the variance due to one update step, ct^(c'). 
To leading order as — > oo, we get 



noU{s,j)[l-UisJ)] 
[1 - 9ix[W 



1 - 9{xi) 
c 



dg{x) 



Af[l-g(xi)]2 dx 



(13) 
(14) 
(15) 
(16) 



Eq. (|16|l gives the increment of the variance of c from one 
update step. To get the total variance of c, we need to 
sum over all updates from ni = to the desired value of 
ni. Provided that there is an upper bound k such that 
dg{x)/dx < K for all x, the total variance of c satisfies 



t(c) < ni 



< 



kN 



(17) 



for xi < 1. (Note that 1/[1 — g{x)] is a nondecreasing 
function because g{x) is nondecreasing.) 



The avalanche is initiated with no* 



"-0*' 



no 



and m = 0. The process ends when no-f-ni = A^ 



and we seek the distribution of no or ni when this hap- 
pens. According to Eq. l(T7|l . the standard deviation of 
c/N scales like I/Vlv, which implies that both hq/N and 
nQ*/N have standard deviations that scale like 1/\/N. 
{xi has zero standard deviation because ni is incre- 
mented by exactly unity on every update step.) Thus 
in the large system limit, the probability of any member 
of the the ensemble of avalanches stopping is negligibly 
small as long as no* /N is finite, and we may treat c as 
exactly conserved as long as this condition holds. 



Using the initial values xi = and no 



A^ 



which determine c, Eq. Q can be rearranged to give 



'0*' 



no 



iV-n[^ 
q Xi T^V- 

'^l-g{Q) 



(18) 



Noting that Uq/N = 1 — xi — nQ*/N, we see that in the 
large A^ limit, the process continues as long as the strict 
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inequality 



xi > [1- g{xi)]- 



lim 

N^oc 



1 - .9(0) 



(19) 



holds, since the inequality implies that Hq* /N remains 
finite. Moreover, in the large N limit it is impossible 
to reach values of xi for which the inequality has the 
opposite sign, because the process stops when Uq* reaches 
zero. 

Note that because of the zero probability of a node 
remaining undamaged when all of its neighbors are dam- 
aged, we have g{l) = 1, which in turn implies that 
Eq. lfTn|l becomes an equality at xi = 1. If Eq. ifT^ 
is satisfied for all xi < 1, the process will be exhaustive 
in the sense that it will not end with a finite value of 
no/N. If, on the other hand, the inequality changes sign 
for xi above some threshold value, then the process will 
terminate when the threshold is reached. If the left hand 
side of Eq. (|19|l forms a tangent line to the right hand 
side of the expression at some value of xi, the process 
will exhibit critical scaling laws. The critical case for ep 
occurs when the when the tangency occurs at xi — 1. 
Examples of these behaviors are presented below and in 
Section HVI 

As an aside, we note that the SP transition is an in- 
stance of criticality at xi — 0. For the above mentioned 
criterion of criticality to hold at xi — 0, the right hand 
side of Eq. H19|) must have the value 1 and the slope —1 
at xi = 0. Thus, the system is critical with respect to SP 
if limjv^oo n^Q*/N = and 



dg{x) 



dx 



1 - .9(0). 



(20) 



Eqs. and (|20|l immediately give a criterion for criti- 
cal percolation on graphs in which every possible directed 
link (including self-inputs) exists with an independent, 
fixed probability, assuming the conventional choice in 
which damage spreads to a given node with probabil- 
ity p from each of its damaged neighbors. In this case we 
have 



5(x)=^P(i^)[l-(l-px)^], (21) 



K=0 



which yields 



dg(x) 



dx 



= pY,p{k)k 



(22) 
(23) 



Piink is the probability that a link exists connecting the 
two randomly selected nodes and p is the probability that 
damage spreads across that link. At the same time, we 
have (K) — punkN. (Recall that K is only the indegree 
of a node, not the total number of links connected to it.) 
Thus Eq. H23|l . which implies that the critical value of 
p is 1/{K), is consistent with the well-known theory of 
Erdos-Renyi graphs. |3^] 

Eq. H23(l applies for any distribution of indegrees so 
long as (K) is well-defined and the source of each input 
is selected at random (so that the outdegrees are Poisson 
distributed). We note that the latter condition is not 
met by random regular graphs (graphs in which all nodes 
have the same outdegree) because the probabilities of 
two nodes getting an output from the same node are 
correlated. 

SP can also be understood by the theory of Galton- 
Watson processes. If limjv^oo "o*/A^ 0, the update 
described by Eqs. is consistent with a Galton- 

Watson processes that has a Poisson out-degree distribu- 
tion with a mean value 



A 



dg{x) 



1 - 5(0) dx 



(24) 



See References 0, |^ [s^l • See Appendix 10 for more 
details on SP in relation to known results. Cases of tan- 
gencies at intermediate values of xi are beyond the scope 
of the present work. 

Returning to the question of the ep transition, it is 
convenient to change variables once again. We define 
XQfi* = 1 — xi and q{xofi*) = 1 — g{xi). In words, q{x) is 
the probability that a randomly selected node will output 
given that each of its inputs has the value with prob- 
ability x. We refer to q as the damage control function as 
it characterizes the probability that damage will be pre- 
vented from spreading to a single node. Equation H19|) is 
then transformed to 



2^0,0* > 9(2^0,0*) 



1 — lim 



9(1) 



(25) 



Critical ep is found when the left hand side of Eq. H25|l 
forms a tangent line to the right hand side of the expres- 
sion at a;o,o* = 0. At criticality, the right hand side of 
Eq. should have the value and the slope 1. Hence, 
the conditions g(0) = and 



dq{x) 



9(1) 



dx 



are required for an ep transition. 



(26) 



This result is closely related to the well-known criterion 
for the presence of a percolating cluster in an Erdos- 
Renyi graph: percolation occurs when the probability Per 
for the presence of a link between two randomly selected 
nodes exceeds 1/N , where N is the number of nodes, [s^ 
In the present context, Per is mapped to punkP, where 



Example: EP on random digraphs 

We now consider the special case of graphs in which 
every possible directed link (including self-inputs) exists 
with an independent, fixed probability. (We have already 
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discussed SP on such graphs.) If damage spreads along 
each directed hnk with probabiUty p, there is no ep tran- 
sition because there is a nonzero probabihty for a node to 
remain undamaged when all of its inputs are damaged. 
A minimal change that allows ep on such graphs is to 
give a special treatment to nodes whose inputs are all 
damaged, in which case the considered node should al- 
ways get damaged. For the same reason, all nodes with 
no inputs must be initially damaged. Other nodes might 
also be initially damaged, and we let this happen with 
a given probability p for each node with at least one in- 
put. For such a network, we can calculate the damage 
propagation function according to 

oo 

g{x) = PiK) [1 - (1 - + (1 - (27) 



K=0 



-{K){l-x)\ 



(28) 



The corresponding damage control function becomes 



q{x) 



1 



(29) 



A necessary condition for the ep transition is derived 
from Eq. H26|l . yielding 



{K)e 



'P(K) 



l-p 



(30) 



For the EP transition to occur, it is also required that 



f{x) = X-q{x){l- p)>Q 



(31) 



for all X e [0, 1] according to Eq. igHJ. If both Eqs. ^ 
and are satisfied, the ep transition occurs at the 
value oip given by Eq. (|^ : 



\tl{K) +ln(l - p) 
W) 



(32) 



Equation H3U|I turns out to be a sufficient and necessary 
condition for the ep transition. Provided that Eq. H3()() 
holds, the first derivative satisfies /'(O) = 0. From the 
observation f"'{x) < 0, it is then straightforward to show 
that f{x) has no local minimum on the interval (0,1). 
Since /(O) = and /(I) > 0, Eq. ^ holds for all x G 

[0,1]. ... 

It is instructive to examine the phase diagram at fixed 
p. A negative value of Pc indicates that the system is 
always in the ep regime, so for (K) < 1 the system 
exhibits ep and it is not possible to observe a transi- 
tion. For (K) > 1, an ep transition can be observed at 
p = Pc- A curious feature of this system is that Pc is not 
a monotonic function of (K) , having a maximum value of 
(1 — p) /e at {K) ~ e/{l — p) and approaching zero as {K) 
approaches infinity. Thus if p is held fixed at any value 
between zero and (1 — p)/e, the system will undergo two 
transitions as (K) is increased from zero. The system 
will begin in the ep regime (i.e. p > pc), undergo a tran- 
sition to subcritical behavior at some (K), then reenter 



0.4 



0.3 



0.2 



0.1 



— 1 1 — I — I — I I I I 



— 1 1 — I — I I I I 



0.0 I ' " ' I ' 




FIG. 3; Phase diagram for ep on random digraphs, where 
damage spreads along each directed link with probability p 
and a node is guaranteed to get damaged in the special case 
that all of its inputs are connected to damaged nodes. All 
nodes with zero inputs are initially damaged, and the other 
nodes are initially damaged with probability p. The gray area 
bounded by a solid line shows the region where ep occurs for 
p = and the dashed lines show the ep transition when p has 
the values 1/4, 1/2, and 3/4, respectively. 



the EP regime for a higher value of {K). The calculated 
phase diagram is shown in Fig. |3 and has been verified 
by direct numerical simulations of avalanches. Roughly 
speaking, at low {K) ep occurs due to the high density 
of initially damaged nodes with no inputs. At high (K), 
on the other hand, ep occurs due to the high probability 
of nodes being damaged because of their large number of 
inputs. 



B. The probability of complete coverage 

An important quantity associated with ep is the prob- 
ability of an avalanche yielding complete coverage of the 
system; i.e., the probability that all sites are damaged by 
the UBA so that u = 0. Let Pcc{N, q;no,no*) denote the 
probability that a UBA on a random network will yield 
complete coverage for a system with a given network size 
N, a given damage control function q, and starting with 
particular values of no and Uq* . For future convenience 
we also define PcciN,q) to be the probability for com- 
plete coverage assuming that each node is initially dam- 
aged with probability 1 — q{l) and we average over the 
corresponding probability distribution for uq*. 

To calculate Pcc{N,q;nQ,nQ*), we note that 



Pcc{N,q;Tn,0) ^0 if m > 0, 



(33) 



since the process stops when uq* = 0. We also have 

Pcc{N,q;0,m) = l for any m, (34) 
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since updating can never create Os. These values of Pec 
can be used for recursive calculation of Pec- Let Uqq* 
denote no + no* , or NxQ f)* . Performing steps OHHl (from 
Section lll A|) one time decreases no o* by 1 as described by 
Eqs. (|5l)-(|Hll. This means that Pcc{N,q;no,nQ*) can be 
calculated for all hq q* = to if Pcc{N, q; no, ng*) is known 
for all hq Q* = TO — 1. The recursion starts at no^o* — 
with Pcc{N, q;0,0) = 1 and uses the boundary condi- 
tions Pcc(iV, g; no,o*, 0) — and Pcc(Af, g; 0, no_o*) = 1 
for no 0* > 0. 

For large N, P^c can be calculated in the framework 
of a continuous approximation. Let p{nQQ*,c) denote 
a continuous version of Pcc{N,q; no, no*). Then, the 
boundary conditions Pcc(-/V, q; J^o.o*! 0) = and Pcc{N, q; 
0, nQ Q*) — 1 are expressed as 



and 



where 



p(%,o*,0) = 1, 



c(2;o,0*) 



"-0,0* 
9(2^0,0* 



(35) 



(36) 



(37) 



In the continuous approximation, the recurrence rela- 
tion that can be derived from Eqs. lO^® is transformed 
to a partial differential equation. In such an update, the 
change no,o* decreases by unity and, for large N, the 
change in c is much less than c itself. In the continu- 
ous approximation, this means that p(no_o*, c) satisfies a 
partial differential equation of the form 



dp 



5no.o* 



hi{no,o*,c) 



dp 
dc 



n'2[nofl*,c)—-^, 



where hi{nQ Q*, c) and /i2('T-o,o*: c) are functions to be de- 
termined. This is recognizable as a ID diffusion equa- 
tion in which no^o* plays the role of time and c the role of 
space. Note that later times in the diffusion equation cor- 
respond to earlier stages of the UBA, since no^o* decreases 
as nodes are converted to Is. The boundary conditions 
on the diffusion are given by Eqs. and l(T7|l . We are 
interested in computing p^Uq q*, c) for values of no o* and 
c corresponding to no* = nf,* and ni = 0. 

The fact that the average of c is constant means that 
the coefficient of the drift term in the diffusion equation 
must vanish; i.e., hi{nQQ*,c) ~ 0. The diffusion coeffi- 
cient, /i2(no_o*j c), is given by 



(39) 



where (T^(c') is the variance of c' when a fixed c is up- 
dated. 

Using Eqs. Hlt)|l and H39|l and converting g's to q's, we 
find 



dp c dq{x) 

driQQ* 2N[q{xQ^Q*)]'^ dx 



The large A'' behavior of Eq. (|40() , with the boundary 
conditions in Eqs. and l|57|) . can be found by ex- 
panding q[x) around x = 0. If q(x) is well-behaved, such 
an expansion can be written as 



q{x) = Oi\x — OLix^ -f 0(2:^). 



(41) 



This expansion can always be performed if the probabil- 
ity P{K) for a node to have K inputs decays as least as 
fast as K~'^ and in the case that px decays slower than 
K~'^ but faster than , only the residue term can be 
affected. See Appendix ^ In particular, the expansion 
is always valid if K has a maximal value. 

The most interesting case in terms of asymptotic be- 
havior is when ai is close to 1 and a2 is positive. With 
suitable TV-dependent transformations of p and its argu- 
ments, described in Appendix^l the large N behavior of 
Eq. 14U|I can be expressed in terms of a function p{t, y) 
determined by the differential equation 



dp 
dt 



1 d^p 



with the boundary conditions 



and 



p{i,l/i)^0 forf<0 



lim p{i, y) = y for y > 0. 



(42) 



(43) 



(44) 



The Crank-Nicholson method can be used to calculate 
p{i,y) numerically in an efficient way. (See, e.g., js^-) 
Appendix El shows that the probability for complete 



(38) coverage is given by 



Pcc(A^, q) « 7V-i/3p[0,7Vi/3(l -ai)], 



(45) 



where N = aiN/a2- The calculation assumes that the 
avalanche is initiated on the nodes whose outputs are 
independent of their inputs, as accounted for in g(l). 

To our knowledge, the critical point for EP has not 
been investigated previously in its own right. Two spe- 
cial cases have been studied, however. First, results for 
numbers of frozen and unfrozen nodes in critical RBNs 
can be mapped to an EP process, as discussed in Sec- 
tion ^3 In this context, frozen nodes in the network are 
considered to be the damaged nodes of the UBA, and the 
scaling with N of the number of unfrozen nodes at the 
phase transition has been investigated for certain class of 
RBNS 25, 29]. 

Second, in the special case that q(x) — x, the exact 
result 



Pcc(iV,a; 1-^ a;; no, no*) 



no* 



no + no* 



(46) 



is obtained. [See Eq. (|B28|) in Appendix IBI] This means 
that the probability for complete coverage is exactly 
n' H./-/V. The simplest realization of q{x) = x is provided 



10 



by a network of one-input nodes with rules that copy the 
input state. Such networks have strong connections to 
random maps from a set of N elements into itself. A 
map T is derived from a network of one-input nodes by 
letting each node map to the node from where its input is 
taken. In this picture, the damage originating from one 
initially damaged node i, corresponds to the set of nodes 
j such that T^{j) = i for some A: > (where denotes 
the fcth iterate of T). Such a j is called a predecessor 
to i. See, e.g., Ref. [3J for an overview of the theory of 
random maps and see Refs. [s^l for results on pre- 
decessors in random maps. See Appendix IeI for analytic 
results that relate UBA to random maps. 



C. On the number of damaged nodes in random 
networks 



In the Sections |lII Al and llll BI we focused on determin- 
ing the parameters that lead to ep (a vanishing fraction of 
undamaged nodes large N limit) and on the probability 
that the number of undamaged nodes will be exactly zero 
(complete coverage). We now consider the full probabil- 
ity distribution for the number of nodes damaged in an 
avalanche in a manner that provides a suitable base for 
understanding both SP and ep in random networks. The 
calculational strategy involves considering a given set of 
n nodes to be the damaged set and computing the prob- 
ability that this is both consistent with all of the Boolean 
rules and the probability that the avalanche will actually 
cover the whole set. The probability of consistency is 
calculated via elementary combinatorics. The probabil- 
ity of reaching the whole set is precisely the probability 
of complete coverage for an avalanche on the sub-network 
of n candidate nodes. For this we can directly apply the 
results of the last section. For the purposes of explaining 
the calculation, we refer to the selected set of n nodes as 
the candidate set. 

We let Pn,N{n) denote the probability that n nodes will 
be damaged in an avalanche, averaged over the ensemble 
of A'^-node networks with a rule distribution character- 
ized by a given damage propagation function g or the 
corresponding damage control function q. We assume 
that the avalanche is initiated by randomly selecting ^ 
nodes to set to 0*, regardless of their Boolean rules, then 
setting to 0* all nodes with rules that always output 1 for 
any inputs. The set of H. initially damaged nodes must be 
a subset of the candidate set. The probability that the 
candidate set contains all of the nodes with "always 1" 
rules will be taken into account by the value of g(0) in 
the expression below for the consistency probability. We 
use the notation (™) for the usual binomial coefficient 
(the number of combinations of k objects chosen from a 
set of m objects). 

The probability P-aM(n) can be expressed as 

Pn,^(n) = {^^l^\p,{n,l-N)PlXn,i-N), (47) 



where P^{n,t,N) and P^^{n,l-N) are defined below. 
The binomial factor counts the number of different sets 
of n — £ nodes that could be damaged in a process cor- 
responding to a given set of £ nodes that are initially 
damaged without regard to their rules. Pc{n, t, N) is the 
probability that a given choice of n — £ nodes assumed 
to be damaged by the avalanche will constitute a final 
state that is consistent with the Boolean rules for each 
node, including the nodes that are initially damaged be- 
cause their rules require it. i'cc("' ^' ^) probabihty 
that the avalanche will not die out before damaging all 
n nodes. This factor is necessary to avoid counting final 
states that contain loops of damaged nodes consistent 
with the rules but unreachable because damage cannot 
spread to the loop from any nodes outside the loop. 

Consistency with the Boolean rules requires that the 
given set of n — £ nodes damaged in the avalanche have 
inputs that cause them to be damaged. In a random net- 
work, the probability that any single node will be dam- 
aged is g{xi), where xi is the fraction of damaged nodes. 
Similarly, the probability that any node will not be dam- 
aged is 1 — g{xi). We are considering candidate sets of 
damaged nodes with xi — n/N . Thus we have 



Pe(n, A^) = [g{n/N)r-'[l ~ g{n/N)] 



(48) 



The computation of Pcc('*i ^! involves the rule dis- 
tribution on the restricted network formed by the candi- 
date set with all inputs from the undamaged nodes re- 
moved. This distribution, g^{x), is different from g{x) 
because Pc already accounts for rules that are not con- 
sistent with the pattern of damage. Thus the spreading 
of damage on the n-node network involves g{nx/N), the 
probability that a rule outputs 1 when a fraction x of the 
n-node candidate set is damaged. The probability must 
be normalized such that it goes to unity when x goes to 
1. (We know that a node in the n-node set should get 
damaged if all of its inputs are damaged.) Thus we have 



9{nx/N) 
gin/N) 



(49) 



or, equivalently. 



1 q[u/N+{l-u/N)x]-q{u/N) 
iN.^u {^) 1 - q{u/N) ■ 

(Recall that u = N — n is the number of undamaged 
nodes after an avalanche.) 

There are two cases of interest for the probability of 
complete coverage of the candidate set. For ep, g{0) > 
and the fixed number £ of nodes arbitrarily selected for 
damage is irrelevant compared to the finite fraction of 
nodes with rules that produce damage for any combina- 
tion of inputs. In this case, we assume £ = 0, which 
allows reduction of Pec to its two-argument form defined 
at the beginning of Section IIII Bl 



Peyn,Q;iV) = Pcc(n,gljv-n)- 



(51) 
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For SP, we have 17(0) = so the avalanche must be 
initiated with a nonzero value of In this case we have 

PlM I- N) = PccXn, g^_^^„; n -£,£). (52) 

Note that P^^{n,t,N) depends on N only through . 

For notational convenience, we now let P^^. stand for 
whichever expression on the right-hand side of Eqs. H51|) 
or H52(l is relevant, and we use u where N — n would be 
the strictly proper form. By combining Eqs. 147|l and 
we get 

PnMn) = [g{n/N)r-'[l - gin/NTP^,. 

(53) 

To make some important features of Eq. (|53|) apparent, 
we introduce the functions 



and 



p{n) 
T{n, k) 

G{x) 



e"n! ' 



n'^{n ~ fc)! ' 



ix)Y fl - g{x) 



X J \ 1 — X 
Then Eq. (|53|l can be rewritten as 

, ,_ pin)p{u) T{n,t) ( n/N \' 
piN) r{N,i)[gin/N)) 

X [G{n/N)rP^^,, 

Stirling's formula, 
yields 



(54) 
(55) 

(56) 



and 



P{n)p{u) _ 1 



(57) 



(58) 



(59) 



(60) 



The factor T{n,£)/T{N,£) is approximately 1 for large 
n and satisfies 



t{N,£) 



< 1 



(61) 



for n < N, with equality if 71 = or ^ = 1 or £ = 0. 
The only factors in Eq. (|57|l that can show exponential 
dependence on TV are the G and P^c factors. Because 
Ppj, is a probability (and therefore cannot exceed unity) 



and G{x) < 1 with equality if and only if g{x) = x, 
Pn,N{n) vanishes exponentially as N goes to infinity for 
any fixed n/N such that g{n/N) ^ n/N. This is con- 
sistent with the above result that the probability of an 
avalanche stopping with xi 7^ g{xi) is vanishingly small. 
[See Eqs. ^ and fl^.] 

For EP, we are interested in the number of undamaged 
nodes, u. We let 



Pn,N{u) = Pn,NiN - u) 



(62) 



and 



Q{x)^G{l-x) 
l-x 



l-x 



q{x) 



(63) 



For EP, 5(0) > and a fixed £ is irrelevant when N ^ 00. 
Hence, we let £ = and rewrite Eq. H57|) as 



Pu,n(u) 



P{n)p{u) 
p{N) 



[Q{u/NrP^,, 



(64) 



To some respects, Pu,7V is similar to Pn,N- the factor 
[p{u)p{n)]/ p{N) is fully symmetric with respect to inter- 
change of n and u; and the role of G{n/N) in Eq. I|57() 
is identical to the role of Q{u/N) in Eq. H64|l . However, 
the behavior of P}.^. for 71 <C iV given by Eq. H52|l is sig- 
nificantly different from the behavior of P},^. for u N 
given by Eq. H51() . 

For EP, we consider damage control functions q{x) that 
can be expanded according to Eq. (|^ . For supercritical 
EP, with ai < 1, Pu,n{u) decays exponentially with u. 
In Appendix ID II we demonstrate that 



lim Pu,A'(u) = (1 - "i) 



N^oo ' • ■ - ' U'. 



/27r 



For critical EP, Eq. gives 

P^,(n,0;iV)=Pcc(n,g^^^_„) 

wn-i/3p[0,ni/3(l-a})]. 



(65) 
(66) 

(67) 



where n = a\n/a2 and a\ and a2 are the first two coeffi- 
cients of the power series expansion of q^{x) about a; = 0. 
With ai — 1 and a2 > 0, a Taylor expansion of logQ(.T) 
about a; = gives 



Q{x) « exp 



v2^3 



(68) 



for small x. This yields that the typical number of un- 
damaged nodes, u, scales like N'^^^. In Appendix ID 21 
we derive the asymptotic distribution of u for large N . 
With u = N~^/^u = {a2/Ny/^u, we find that the large 
N limit of the probability density for u is 



Piu) 



exp(^ 



p(0,2u). 



(69) 
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Eq. (|57() is suitable for understanding SP as well as 
EP. For SP, g{0) = and £ > 0. In the large TV limit, 
SP is a branching process with a Poisson distribution in 
the number of branches from each node. The average 
number of branches per node is given by the derivative 
of g{x) at X = 0, because limx~>o g{x)/x is the average 
number of nodes that will be damaged in one update as 
a direct consequence of damaging a single node in the 
large network limit. In Appendix [nj we re-derive known 
results on SP in the framework of our formalism. 



IV. AN APPLICATION: FROZEN NODES IN 
RANDOM BOOLEAN NETWORKS 

An important application of our results on ep in ran- 
dom networks is the determination of the size distribu- 
tion for the set of unfrozen nodes in 2-input random 
Boolean networks, a subject of interest since the intro- 
duction of the Kauffman model in 1969 j^. The Kauff- 
man model was originally proposed as a vehicle for study- 
ing aspects of the complex dynamics of transcriptional 
networks within cells. 

In a Boolean network, there are usually some nodes 
that will reach a fixed final state after a transient time 
regardless of the initial state of the network. For most 
random Boolean networks, nearly all of these nodes can 
be found by a procedure introduced in Ref. and ap- 
plied numerically in Ref. [S^l • We refer to nodes identified 
by this procedure as frozen. 

The nodes that cannot be identified as frozen are la- 
beled unfrozen. Their output may switch on and off for 
all time or simply have different values on different at- 
tractors of the network dynamics. A frozen node will 
always reach its fixed final state regardless of the ini- 
tial state of the network. The converse is not true: an 
unfrozen node can have a fixed final state that is indepen- 
dent of the initial state due to correlations that are not 
accounted for in the identification procedure for frozen 
nodes. In a typical random Boolean network, the num- 
ber of nodes that are mislabeled in this sense is negligible 
|27j . For the purposes of investigating dynamics of the 
network at long times, one is interested in the size of the 
unfrozen set. 

The procedure for identification of the frozen nodes 
starts by marking all nodes with a constant output func- 
tion as frozen. There may then be nodes that, as a conse- 
quence of receiving one or many inputs from frozen nodes, 
will also produce a constant output. These nodes are also 
marked as frozen, and the process continues iteratively 
until there are no further nodes that can be identified as 
frozen. 

We note here that the process of finding frozen nodes 
in a RBN can often be framed as a UBA, where the prop- 
erty of being frozen corresponds to damage. That is, the 
process of identifying frozen nodes involves continually 
checking all nodes to see whether their inputs are frozen 
in such a way that they themselves become frozen, a pro- 



cess which satisfies the conditions for UBA. The dam- 
age propagation and damage control functions for the 
UBA are determined by the relative weights of different 
Boolean logic functions in the rbn. By changing these 
weights, one can observe a transition in the dynamical 
behavior of rbns corresponding precisely to the ep tran- 
sition in the UBA. We consider here rbns with exactly 
two inputs at each node, with some explicit choices for 
the weights of the Boolean logic functions that permit 
observation of both sides of the transition. 

The only restriction required for mapping the freez- 
ing of nodes in a rbn to a UBA system is that the logic 
functions in the RBN be symmetric with respect to the 
probability of freezing being due to true and false in- 
puts. That is, the probability that a node with a certain 
set of frozen inputs will itself be frozen should not de- 
pend on the values of the frozen inputs. This condition 
is satisfied for the most commonly investigated classes 
of rule distributions, where there is a given probability 
p for obtaining a 1 at each entry in the truth table for 
each rule. If the above mentioned symmetry condition 
were violated, it would be necessary to distinguish nodes 
frozen true from nodes frozen false, which would mean 
that the state of a node could not be specified by a bi- 
nary variable. For the rest of this section we consider 
only RBNS that respect the symmetry condition. 

It is useful to distinguish different types of Boolean 
logic functions. A canalizing rule is one for which the 
output is independent of one of the inputs for at least 
one value of the other input. Among the 16 possible 2- 
input Boolean rules, 2 rules are constant ("always on" or 
"always off"), 12 rules are non-constant and canalizing, 
and 2 rules are non-canalizing (xOR and not-xOR). The 
original version of the Kauffman model assumes that all 
2-input Boolean rules are equally likely, which turns out 
to give critical dynamics. 

Let Pi denote the probability that a randomly selected 
node's output is frozen if exactly i of its inputs are frozen. 
The damage propagation function g(x) can be expressed 
directly in terms of pf. 

g{x) = po{l - xf + 2pix{l ~x)+ p2x'^. (70) 

Nodes with constant rules are guaranteed to be frozen. 
(These nodes will initiate the UBA.) Nodes with non- 
constant canalizing rules are unfrozen if both inputs are 
unfrozen, and they are frozen with probability 1/2 if ex- 
actly one randomly selected input is frozen. Nodes with 
rules that are non-canalizing become frozen if and only 
if both of their inputs are frozen. Finally, if both inputs 
are frozen, the output of any 2-input rule is frozen. Thus 
for the 2-input Kauffman model, po — 1/8, Pi — 1/2, and 
P2 = 1- 

If the two non-canalizing rules in the 2-input Kauffman 
model are replaced by canalizing rules, pi becomes 9/16, 
whereas po £^nd p2 are unchanged. Such networks exhibit 
supercritical ep. To get a subcritical network, we replace 
two of the canalizing rules with non-canalizing rules and 
get pi = 7/16. (Note that some care must be taken to 
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FIG. 4: The functions (a) g{x) = 1 - q{l - x) and (b) 
G{x) = Q(l — x) for three 2-inputs rule distributions. All 
three distributions have po — 1/8 and p2 = 1, whereas pi 
takes the values 7/16, 1/2, and 9/16. The case that has 
pi = 1/2 (marked with =) is critical with respect to ep and 
corresponds to the propagation of frozen node values in the 
original Kauffman model. The other cases pi — 7/16 (<) and 
Pi — 9/16 (>) are subcritical and supercritical, respectively. 
The dashed line in (a) shows the identity function a; t-» s. 



maintain the true-false symmetry mentioned above.) 
The functions g{x) and G{x) for critical, supercritical, 
and subcritical rule distributions are shown in Fig. 01 

As can be seen from Fig.01 a small change in g{x) may 
lead to a qualitative change in G{x) for rule distributions 
close to criticality. Such changes have a strong impact 
on the avalanche size distribution for large N. Figure [S] 
shows the probability density distribution of the fraction, 
n/N, of nodes that are affected by avalanches in networks 
with the above mentioned rule distributions. The proba- 
bility distributions are obtained by recursive calculation 
of the distribution of no* as ni increases. The recurrence 
relations are obtained from Eqs. ®-(|HJl and the result 
is exact up to truncation errors. To verify these calcula- 



FIG. 5: The probability density distribution NPn,N{n) 
with respect to the fraction of nodes {n/N) involved in an 
avalanche. The rule distributions have the same g[x) as dis- 
played in Fig. 13] showing rule distributions that are (<) sub- 
critical, (=) critical, and (>) supercritical with respect to 
EP. The displayed networks sizes, A^, are 10 (large dots), 100 
(small dots), 10^ (bold line), lO", 10^ and lO'' (gradually 
thinner lines). 



tions, we generated 10^ random Boolean networks of size 
N — 10"^ for each of the above described rule distribu- 
tions. The distributions in the numbers of frozen nodes 
in those networks are displayed in Fig. [S] 

In Fig. [71 the probability distributions of the number 
of undamaged nodes, u, are shown in comparison to the 
asymptotic results in Eqs. H65() and H69() . Our analytic 
results are strengthened by the data in Fig.[7|as the dis- 
tributions for finite networks approaches the predicted 
asymptotes. Finite size effects are clearly visible in the 
critical case even for network sizes as big as TV = 10^, 
whereas convergence in the supercritical case is achieved 
for iV > 10^. 

Kaufman, Mihaljev, and Drossel studied distributions 
of unfrozen nodes in 2-input critical RBNs using a method 
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FIG. 6: A numeric comparison between analytic calculations 
(black lines) and explicit reductions of random Boolean net- 
works (gray lines). For both cases, the probability density 
distribution NPn,N{ii) is displayed as a function of n/N. The 
rule distributions have the same g{x) as displayed in Figs. |2 
and|3 showing rule distributions that are (<) subcritical, (=) 
critical, and (>) supercritical with respect to ep. The UBA 
rule distributions are realized in random Boolean networks 
by rule distributions with the following respective selection 
probabilities: 1/8, 1/4, 5/8 —pr, and pr for a constant rule, 
a rule that depends on exactly 1 input, a canalizing rule that 
depends on 2 inputs, a 2-input reversible rule. The values of 
Pr axe (<) 0, (=) 1/8, and (>) 1/4. For each rule distribution, 
10^ networks were tested. 



similar to ours in that differential equations for popula- 
tions of different types of nodes are developed from a 
discrete process in which frozen nodes are identified by 
the propagation of information from their inputs [29l |. 
Their result for the numbers of unfrozen nodes in 2-input 
critical RBNs corres pon ds to a particular application of 
Eq. (|69|l . In Ref. [23, the function corresponding to 
P{u) [which they call G{y)] is determined by running a 
stochastic process and a numerically motivated approxi- 
mation is proposed: 




Probability 




P{u) « 0.25 exp(- 



l-0.5\/a + 3w 



(71) 



200 



FIG. 7: Rescaled versions of the probability distributions 
displayed in Fig. |3 (a) the probability density for the criti- 
cal case, with respect to the rescaled number of undamaged 
nodes, u = {a2/Nf/^u = u/{4:N^^^); (b) the probability dis- 
tribution Pu,n{u) for the supercritical case without rescal- 
ing. The displayed networks sizes. A*', are 10 (large dots), 
100 (small dots), 10^ (bold line), lO", 10^ and 10^ (gradu- 
ally thinner lines). The analytically derived asymptotes are 
shown as dashed lines. In (b), the distributions for networks 
of sizes 10", 10^ and 10'^ are not plotted because they are 
indistinguishable from the asymptotic curve. 



The scaling law P{u) cx u^^/^ for small u is also derived 
analytically in Ref. 29]. 

For large x, Eqs. (|42|1 - H44() imply p{0, x) (x x for large 
positive X. This means that 



P(i2) « y ^ exp(- 



(72) 



for large u. Thus the large u limit of Eq. (|71|l differs from 
the exact result by a factor of (3/4)-\/7i'/2, an underesti- 
mate of about 6%. 

We are able to improve further on Eq. (|71|l by nu- 
merical investigations of p(0, x) calculated by the Crank- 



Nicholson method (see, e.g., f^) using Eqs. 

We find that the high-precision numerical results are fit 

by the function 



Piu) 



exp(- 



) 1 + 



1 



3.248m- 



4.27^2 4.76u3 
(73) 



with a relative error that is maximally 0.25% and van- 
ishing for large u. 

By explicitly keeping track of the populations of nodes 
with each of the different types of Boolean logic functions 
as links from frozen nodes are deleted, Kaufman, Mihal- 
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jev, and Drossel |29j also derive results for other quan- 
tities, such as the number of links in the sub-network 
of unfrozen nodes. The ep formalism described above 
can be applied once again to investigate these additional 
quantities in a broader class of networks. Detailed re- 
sults for RBNs with various degree distributions will be 
presented elsewhere. 

V. SUMMARY AND DISCUSSION 

Unordered binary avalanches can in some cases lead 
to damage on every node or almost every node of a net- 
work, a phenomenon we have dubbed exhaustive percola- 
tion. We have studied a broad class of random networks 
that can exhibit ep. We have shown how to calculate 
the probability Pec (N) that complete coverage occurs (i.e 
that all nodes are damaged) and also derived expressions 
for the probability distribution P{u) of the number of 
undamaged nodes, u, in the large N limit when ep does 
occur. A logical curiosity in our approach is the fact 
that the calculation of P{u) involves application of the 
Pec result to subnetworks containing candidate sets of 
damaged nodes. 

Our primary results flow from the realization that all 
of the relevant information about a UBA defined on a 
random network is contained in the damage propagation 
function g{x) or, equivalently, the damage control func- 
tion q{x). We derive scaling law exponents and exact 
results for the distribution of u that are valid for a broad 
class of random networks and Boolean rule distributions 
in the ep regime and for networks at the ep critical point. 
This class includes the UBAs that determine the set of 
frozen nodes in RBNs with more than two inputs per node 
and therefore constitute a generalization of the results on 
the set of unfrozen nodes in rbns presented in Ref. 1291 . 
Interestingly, the asymptotic behavior found in Ref. |29| 



for the distribution of u at the critical point is shown to 
be valid for a broad class of network problems. 

For networks outside the above mentioned class but 
within the framework of UBA, we find connections to pre- 
vious work on Galton- Watson processes and random 
maps [3^ . The central result of our investigations is dis- 
played in Eqs. Ht)t)|l and Ht)9|) . which provide explicit for- 
mulas for the probability of finding u undamaged nodes 
after an avalanche runs to completion. The out-degree 
distributions of the networks described by our formulas 
are all Poissonian, but the in-degree distributions may 
have different forms, including power laws, so long as 
the probability of having in-degree K decays faster than 
K^^. The exact nature of the EP transition on networks 
with broader in-degree distributions is an interesting is- 
sue for future research. Further work is also needed 
to handle correlations between input links to different 
nodes, a situation that arises, for example, in random 
regular graphs or networks with scale free out-degree dis- 
tributions. 

Our original motivation for studying EP arose from at- 
tempts to understand the dynamical behavior of rbns. 
We have described one nontrivial example of how the ep 
formalism is relevant: the calculation of the probability 
distribution for the number of unfrozen nodes in any RBN 
with a rule distribution that leads to a given damage con- 
trol function q for the associated UBA. The problem of 
determining how many of the unfrozen nodes are actu- 
ally relevant for determining the attractor structure of 
the RBN can also be framed as an ep problem, which will 
be addressed in a separate publication. 
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APPENDIX A: CALCULATION OF THE 
DAMAGE CONTROL FUNCTION 

Let pk denote the probability that a rule has K inputs, 
and let Pq^K^ m) denote the probability that the output 
value is zero of a rule with K inputs fed with m zeros 
and K — m ones. Then, the the damage control function 
is 



K=0 m=0 



q{x)^ 5]pK^Po(i^,m)(^)x"(l-x)^ 



Eq. IjAlfl can be written as 

q[x) = ao + aix + a^x^ + 

where 



(Al) 



(A2) 



K=i m=0 



K 



K — m 
i ~ m 



(A3) 



The expansion in Eq. I|A2|I is well-defined up to the 
first term such that the sum in Eq. (|A3() is not abso- 
lute convergent. The factor scales like for 
large K and Po{K,m) < 1. Hence, Ui is well-defined if 
Sk=o ^^Pk is convergent and this is true if pk decays 
faster than K~^~'^ . 



In addition, the requirement that the output of each 
rule in the rule distribution is 1 if all of its inputs have 
the value 1, yields that oq = 0. Thus, the expansion 



q{x) = aix — a2X^ + 0{x^), 



(A4) 



is valid for all rule distributions such that pk decays 
faster than K^"^. In the case that pK decays slower than 
K~* but faster than K~^, only the residue term can be 
affected. 



APPENDIX B: PROBABILITY FOR COMPLETE 
COVERAGE 

Here, we assume that the expansion in Eq. ljA4|l is 
well-defined. Then, we get 



dp 



d^p 



dn, 



0,0* 



[1 + O{no,o*/N)], (Bl) 



and 



.{x)/N 



a2 



(B2) 



To remove the dependence of riQ.o* from the leading 
order term of the diffusion rate in Eq. (jBip . we let t = 
— l/no,o*- By also letting y = 1 — aic/N, we rewrite 
Eq. ljBl|) to a form that easily can be rescaled as N grows. 
We get 



dp 
'dt 



1 ~ y d'^p 
2 



and 



Vmin 



a2 
aiNt 



(B3) 



(B4) 



where y = j/min is the transformed value of Cmax- The 
boundary conditions are p = Q ior y — j/min and p — 1 
for y ~ 1. 

The N dependence of the leading order term of the 
boundary condition in Eq. ljB4|) can be removed by rescal- 
ing of y and t. Typically, a2 > and this is the case that 
we will focus on. [Note that Eq. H25|l means that a2 
must be nonnegative at the transition.] If a2 — 0, either 
q{x) = a; or q{x) = x — Umx"^ + • ■ ■ with m > 2 (apart 
from some pathological special cases). The first case, 
q{x) = X, is a special case that is convenient for analytic 
calculation, whereas the latter case require calculations 
analogous to the calculations for a2 > 0. We will come 
back to the case q{x) = x. 

For Q!2 > 0, we rescale y and t according to 



5 = 



y 



V 



and 



t = N'^/^t, 



(B5) 



(B6) 
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where 



N=^N. 



(B7) 



Then, 



where 

ymin = f-^[l + 0(^-)]- 



(B9) 



The boundary conditions arc p = for ?/min and p = 1 
for y = N^l"^ . The only plausible limit of p as i ^ — oo is 
•p = y/N^/^. To get a motivation that is mathematically 
acceptable, one needs to relate the original integer based 
formulation of the problem in Eqs. (0)-®- The large N 
behavior of Eq. ||SJ|, 



yields 



lim U{sJ) = 1/^0,0*, 



lim PcciN,q;no,no*) ^ 



(BIO) 



(Bll) 



Eq. (|Blip can be shown via induction. The induction is 
initiated by 



lim Pcc(iv,<z;i,0) = 



and 



N 



lim P,:,c{N,q;0,l) = 1, 

N^oo 



(B12) 



(B13) 



which means that Eq. (|B11|) is true for Uq q* = 1. To 
obtain the induction step, we assume that Eq. (|Blip is 
true for = uq q* — 1. Then, we get 



lim Pcc{N,q;n'o,n'Q*) ^ 



which leads to 



lim Pcc{N,q;no,no*) = 



no* + ?^o/no,o* - 1 



710,0* - 1 



Hqq* 



(B14) 

(B15) 

(B16) 
(B17) 



that completes the induction step. Eq. I|B11|I means that 
the value of p approaches a linear function of y for t ~ 
N^^^/uq as N oo. Hence, the boundary condition for 
t — oo is p = y/N^/^. 
Rescaling of p according to 



(B18) 



gives the boundary condition 

lim p — y. 



(B19) 



If Eq. (|B19(I is extended to be valid for all non-negative 
y, the boundary condition at y = iV^/^ can be dropped. 
In the limit of large N, Eq. (|B8|I becomes 

dp 1 d'^p 
dt " 20^' 



(B20) 



With p is written on the form p{t, y), the boundary con- 
ditions are 



and 



p{i, = 



. lim Pit, y) = y 



for t < 



for y > 0, 



(B21) 



(B22) 



as — > oo. 

The solution to Eqs. HB20p ~ (|B22p can be calculated 
numerically. By expressing the transformed variables t 
and y in terms of more fundamental quantities, we get 



PcciN,q;na,no*) 



N~^/^p 



'^0,0* ' 



Nq{no,o*/N) 



and 



PcciN,q;N ~ n„*,nQ*) 



N-'/'p 



N 



1/3 



1 - 



Nq{l) 



(B23) 



(B24) 



for large N. 

If the avalanche is initiated by letting each node start 
from with probability q{l), we get 



K*) = - 9(1)] 



and 



(7(4*) = V^g(l)[l-g(l)]. 



(B25) 



(B26) 



Provided that q{x) does not depend on iV, N is fixed 
and the spread in y that correspond to the initial value of 
m will go to zero as iV ^ oo. Also, N'^/'^ /N approaches 
zero as A'^ — *■ oo. In this case, the probability for an 
avalanche to yield complete coverage is given by 

Pcc(iV, q) « N-^/^p[Q, N^/\l - ai)], (B27) 

for large N . Only the first two arguments to Pec are kept 
in Eq. (|B27|1 . because the process is fully determined by 
N and q. 

In the special case that q{x) = x for all x G [0,1], 
Eq. (jSJ yields a = l/no,o*j which is a strong form of 
Eq. HBlOfl . By using the same induction steps that lead 
from Eq. IBlOp to Eq. (|B11|) . we conclude that 



Pcc(A^,a; 1-^ a;;no,no*) 



(B28) 
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APPENDIX C: ASYMPTOTES FOR SPARSE 
PERCOLATION 

Provided that the derivative of g(x) is well defined at 
a; = 0, we let A = g'{0), where g'{x) denotes the deriva- 
tive of g{x). Then, 



iV— >oo 



which means that 



(Cl) 



(C2) 



(C3) 



according to Eq. (|B28|I . Thus, the large N limit of 
Eq. ISZI) is 



and 



lim QN^N-nix) = X 



lim P^,,{n,£;N) = - 



Thus, there are two competing processes as N ^ oo. 
The sensitivity with respect to q is limited by the vari- 
ance of the number of nodes with initial state 0*, because 
this variance can be seen as a rescaling of q. The change 
in q{x) by such a rescaling scales like q{x)/^fN for large 
N . If q{x) has a well-defined nonzero derivative at a; = 0, 
the difference q\j ^{x) — q{x) scales like q{x)/N for large 
N. Hence, the decrease in the difference between qj^^ 
and q dominates over the increase in sensitivity, meaning 
that 



Thus, 



lim ^---(^-"'Q'^) = 1. 

AT^oo Pcc{N) 



lini :^i4S = p(7.)A«e"(i-^) 



(D2) 



(D3) 



(D4) 



where (uA)" should be interpreted with the convention 
that O'^ = 1 in order to handle the case u — properly. 



lim PnwH =/9(n)r(n,^)A"-^e"(i-^)- 
N—>oo ' n 

£(nA)"-^_„, 



n{n-£)l 
For large N, Eq. (|C5p yields 

i 



N 



lim PnM{n) 



2tt 



(C4) 
(C5) 

(C6) 



Due to the correspondence to well investigated branch- 
ing processes, Eq. (IC5|) is not a new result. For the special 
case of ^ = 1, Eq. IjCSp is given explicitly in Ref. |33|, and 
the general form of Eq. (|C5|I can easily be obtained by 
the theorem presented in Ref. [s^ . 



APPENDIX D: ASYMPTOTES FOR 
EXHAUSTIVE PERCOLATION 

In analogy with our investigation of SP, we assume 
that q{x) has a well-defined derivative at x = and let 
A = ^'(O). For EP to be likely in the large N limit, it 
is required that q{x) < x for all x, meaning that A < 1. 
The large N behavior of Eq. (|64|l is partly explained by 



lim 



piu)p{N ~ u) 
p{N) 



[Q{u/N)f - p{u)X 



> ti u(l-A) 



(Dl) 



but it remains to investigate the role of Pcd^ ^ 
u,0;N). To this end, we consider the ratio PcdN — 
u,0; N)/ Pcc{N). [Here, we have dropped the argument 
q bom PcciN, q).] 

When N ^ oo, there are two processes that influence 
on this ratio: „ approaches q and N increases. The 
increase of TV makes the involved probabilities more sen- 
sitive for the shrinking differences between qj^ ^ and q. 



1. Limit distributions for supercritical EP 

If < A < 1 and x = is the only solution to q{x) = 
X in the interval < a; < 1, the exponential decay of 
[Q{u/N)]^ , in Eq. with increasing u ensures that 



lim Pu.n{u) = 1 



u=0 



and 



lim PcciN)] y (!^e-"^ = 1. 



(D5) 



(D6) 



u=0 



Thus, limAT^+oo has a unique value for each A. This value 
can be calculated by considering the simplest case, qix) = 
Xx. From the definition of the spreading process, we get 

PcciN, X ^ Ax; no, no*) = Pcc{N,x ^ a;; no, no*). 

(D7) 

Then, Eq. (|B28p and averaging over initial configurations 
yield 

PcciN,x^ Xx) ^ 1- X, (D8) 
which means that 



lim PcciN, q) = 1 - A 



(D9) 



for all q that satisfy the above mentioned criteria. We 

get, 



lim P„_^(u) = (l-A)^^e-"\ 



which for large u means that 
1- A 



lim Pn,Niu) 



"(l-A);^«y-l/2^ 



(DIO) 



(Dll) 
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2. Scaling at the EP transition 

This section aims to derive the asymptotic distribution 
of u for large N for critical EP with ^2 > and I = Q. 
Define aj, a\ and analogous to the definitions of ai, 
a2 and N in Eqs. (HH) and (|B7|) . The derivatives of q^{x) 
in Eq. (|50|l at a; = are given by 



Hence, 



{q')'{0)^q'{ulN) 



l~u/N 
1 - q{u/N) 



(D12) 



and 



and we get 

N^=N + 0{u/N) (D14) 
and 

= ai + (ai - 1 - 2a2)M/A^ + 0{u^/N'^). (D15) 

For critical networks, with ai = 1, we get N — N/a2 
and 



a\ = l- 2u/N + 0{u^/N^). 



(D16) 

Insertion into Eq. (jB27|) yields 

PjjAT - M, 0; N) « iV-i/3p(0, 27V-2/3u) (di7) 
for N. 

Because Q{x) < 1 with equality if and only if q{x) = x, 
Pu.N that vanishes exponentially, as N goes to infinity, 
for any fixed u/N such that q{u/N) ^ u/N . For a typical 
network with 0:2 > at the EP transition, the only solu- 
tion to q{x) = a; is X = 0. For such a network, the large 
TV behavior of Pu.iv is found by expanding Q{x) around 
X = 0. To the leading non- vanishing order, we get 



Q(x) « exp 



,,2^3 



which yields 



Pu.n{u) ~ p(u)exp 



2iV2 



(D18) 



(D19) 



Pu,^(u) « iV-i/3p(u) exp(^-i^p(0, 2N~^^^u), 

(D20) 

with asymptotic equality for large N. The probability 
density, P{u), for the distribution of u = N~'^/'^u as N —> 
00 approaches 



'2ttu 



(D21) 



APPENDIX E: EXACT RESULTS 

A network with g{x) = a; is critical for all x. For such 
a network, G{x) — 1 for all x and P^f^{n,£; N) — £/n. 
Hence, Eq. H57|) yields 



Pn,N{n) 



p{n)p{N - n) T{n,£) £ 
p{N) t{N, £)n 
£ fN -£\ n'^-^iN ^n)^' 
n\n ^ £ , 



For n and N satisfying n':^ 1 and N — n':^ l,we get 

£ Vn 



(El) 
(E2) 



Pn,N{n) 



' Zn n y/n{N - n) 



(E3) 



In the special case £ — I, Pn,A'(") is the distribution 
of the number of predecessors to an element in a random 
map. This distribution, which is consistent with eq. (|E2|) 
for £ = 1, was obtained in Ref. ^.nd restated in Ref. 

For completeness, we provide an explicit expression for 
the distribution of avalanche sizes in the case that g{x) 
is a first order polynomial. From Eq. |S7J), we get 



Pn,N{n) ^Pum{u) 

_ fN - £ 

\ n 



3(0) 



nq{0) + uq{l] 



(E4) 



